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Abstract 

Background: MicroRNAs (miRNAs) are small, endogenous RNAs that play important regulatory roles in 
development and stress response in plants by negatively affecting gene expression post-transcriptionally. 
Identification of miRNAs at the global genome-level by high-throughout sequencing is essential to functionally 
characterize miRNAs in plants. Drought is one of the common environmental stresses limiting plant growth and 
development. To understand the role of miRNAs in response of plants to drought stress, drought-responsive 
miRNAs were identified by high-throughput sequencing in a legume model plant, Medicago truncatula. 

Results: Two hundreds eighty three and 293 known miRNAs were identified from the control and drought stress 
libraries, respectively. In addition, 238 potential candidate miRNAs were identified, and among them 14 new 
miRNAs and 15 new members of known miRNA families whose complementary miRNA*s were also detected. Both 
high-throughput sequencing and RT-qPCR confirmed that 22 members of 4 miRNA families were up-regulated and 
10 members of 6 miRNA families were down-regulated in response to drought stress. Among the 29 new miRNAs/ 
new members of known miRNA families, 8 miRNAs were responsive to drought stress with both 4 miRNAs being 
up- and down-regulated, respectively. The known and predicted targets of the drought-responsive miRNAs were 
found to be involved in diverse cellular processes in plants, including development, transcription, protein 
degradation, detoxification, nutrient status and cross adaptation. 

Conclusions: We identified 32 known members of 10 miRNA families and 8 new miRNAs/new members of known 
miRNA families that were responsive to drought stress by high-throughput sequencing of small RNAs from M. 
truncatula. These findings are of importance for our understanding of the roles played by miRNAs in response of 
plants to abiotic stress in general and drought stress in particular. 




Genomics 



Background 

Endogenous small interfering RNAs (siRNAs) and 
microRNAs (miRNAs) are the two most abundant 
classes of plant small RNAs (sRNAs). The small RNAs 
are processed in the nucleus from longer precursor tran- 
scripts that form distinct secondary structures. The 
miRNAs down-regulate gene expression by targeting 
specific messenger RNAs (mRNAs) in both plants and 
animals [1,2]. 

miRNAs were initially discovered in Caenorhabditis 
elegans as developmental timing regulators in 1994 [3]. 
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The existence of miRNAs in organisms including plants, 
mammals and virus has widely been recognized. The 
biogenesis of miRNAs in plants is a multi-step enzy- 
matic process involving incorporation of miRNAs into 
the RNA-induced silencing complex (RISC), and then 
miRNAs bind target mRNAs to direct for cleavage miR- 
NAs with near perfect complementarity and/or inhibit- 
ing translation of those with lower complementarity 
[4-7]. There has been ample evidence demonstrating 
that miRNAs play a regulatory role in diverse biochem- 
ical and physiological processes in plants [1,4]. For 
instance, miRNAs have been shown to play a role in the 
modulation of the processes associated with growth and 
development in plants, including leaf morphogenesis, 
floral organ and root development [6,8-13]. In addition. 
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recent studies also revealed that miRNAs are involved in 
responses of plants to various abiotic and biotic stresses. 
These include drought [14-18], cold [19-21], salinity 
[22], nutrient starvation [23-26], oxidative stress [27], 
submergence [28], UV-B radiation [29,30] and virus 
[31,32]. 

Plants are hardly grown under optimal conditions, and 
have to frequently suffer from drought stress due to short- 
age of water supply. Not unexpectedly, the involvement of 
miRNAs in response of plants to drought has been evalu- 
ated in several plant species. A number of miRNAs have 
been identified to be associated with response to drought 
stress in several species such as Arabidopsis thaliana 
[14,33] and Oryza sativa [17,18]. For instance, Zhou et al. 
(2010) identified 11 down-regulated miRNAs and 8 up- 
regulated miRNAs in response of rice to drought [18]. It 
has been reported that miR169g is the only member 
induced by drought in the miR169 family of rice and its 
expression is regulated by drought [17]. Li et al. (2008) 
found that miR169a and miR169c are substantially down- 
regulated by drought, leading to the enhanced resistance 
to drought in Arabidopsis because one of the miR169's 
targets, NFYA5 (Nuclear Factor YA5), is a crucial tran- 
scription factor regulating the expression of a number of 
drought stress-responsive genes [14]. Recent studies also 
demonstrated that miRNAs play a role in control of 
drought resistance in maize by modulating the expression 
of MAPK (mitogen- activated protein kinase), PLD (phos- 
pholipase D), PDH (praline dehydrogenase) and POD 
(peroxidase) [16]. More recently, Trinadade et al. (2010) 
reported that miR398 and miR408 are up-regulated by 
water deficit in Medicago truncatula, leading to a down- 
regulation of their target genes of COXSb, CSDl and plan- 
tacyanin [15]. In Triticum dicoccoides, 13 miRNAs are 
differentially regulated in response to short-term drought 
stress by miRNA microarray chip [34]. Despite of numer- 
ous studies on the involvement of miRNAs in abiotic 
stress in general and drought stress in particular, there has 
been no report on the systemic identification of drought- 
responsive miRNAs in leguminous plants at the global 
genome-level by high-throughout sequencing. 

Leguminous plants account for one third of primary 
crop production in the world and are important sources 
for the consumption of human and animals [35]. Drought 
stress is one of the most frequently occurred environmen- 
tal stresses that limit crop yield world-wide. Medicago 
truncatula is an annual legume species distinguished by 
its small diploid genome and easy transformation, and is a 
model plant to study functional genomics of legume plants 
[36] . Although numerous miRNAs have been identified in 
various plant species and the mechanisms underlying their 
actions are being unravelled, the discovery of new miRNAs 
in plants on a genome-wide scale is essential for fianctional 



characterization of miRNAs. The traditional sequencing of 
relatively small-size cDNA libraries of plant sRNAs from 
Arabidopsis, rice and poplar with Sanger method has led 
to the conclusions that plant miRNAs are highly con- 
served [37] . However, a small number of miRNAs was not 
detected in the genomes of some species, indicating that 
miRNAs may evolve more recently [38]. In addition, the 
non-conserved miRNAs are often expressed at a lower 
level than the conserved miRNAs, thus many non- 
conserved miRNAs cannot be detected in small-scale 
sequencing studies. In this context, high-throughput 
sequencing has been used to identify non-conserved miR- 
NAs in several species [20,39-53]. To understand the role 
of miRNAs in response of plants to drought stress, we 
identified a number of conserved and non-conserved miR- 
NAs that were responsive to drought by high-throughput 
sequencing, and their potential role in drought tolerance 
was discussed. 

Results 

High-throughput sequencing of small RNAs from 
M. truncatula 

Two sRNA libraries in which M. truncatula seedlings were 
treated with drought stress (DS) and without stress, control 
(CK) constructed from shoots of seedlings grown under 
the two treatment-regimes were sequenced by an lUumina- 
Solexa sequencer. High-throughput sequencing of CK and 
DS libraries led to generation of 6,808,877 and 6,874,742 
primary reads, respectively. There were 6,420,234 clean 
reads (2,685,754 unique sequences) for CK and 6,505,965 
clean reads (2,937,417 unique sequences) for DS after 
initial processing (Table 1). The length distribution of 
reads showed that the majority of the reads was 20-25 nt 
in size, of which the class of 24 nt was the most abundant 
group followed by the group of 21 nt class (Figure la). 

The common/specific sequences were analyzed 
between CK and DS libraries for the total and unique 
sequences. There were 19.14% and 44.95% specific 
sequences in the DS library for total and unique 
sequences, respectively. 

After initial processing, the high-quality small RNA 
reads were mapped to the M. truncatula genome 
sequence (Mt3.0) using SOAP [54]. The number of total 
and unique sequences that matched genome was 
4,875,034/1,713,480 and 4,819,359/1,867,375 in the two 
libraries, respectively (Table 1). 

Identification of known miRNAs 

The known miRNAs in M. truncatula were identified by 
comparing with the up-dated miRNAs database (miR- 
Base 17, released in April 26, 2011) [55]. We identified 
283 and 293 known miRNAs in the two libraries, 
respectively (Table 1). Sixty-six members belonging to 
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Table 1 Statistics of sRNA sequences for control (CK) and drought stress (DS) libraries 



Number of reads 



Number of unique sequences 



Control (CK) 

Primary reads 
Clean reads 

Sequences mapped to tlie genome 
matcin known miRNAs 



6,808,877 
6,420,234 
4,875,034 
509,981 



2,685,754 
1,713,480 
283 



Drought stress (DS) 
Primary reads 
Clean reads 

Sequences mapped to tlie genome 
Matcli l<nown miRNAs 



6,874742 
6,505,965 
4,81 9,359 
588,990 



2,937,417 
1,867,375 
293 



14 conserved miRNA families were obtained after 
removing those miRNAs whose expression levels are too 
low to be analysed for differential expression. Of these 
families, the most abundant two reads were miR156 and 
miR166. Among the identified miRNAs, miR399 con- 
tained the most members, including miR399a-i and 
miR399k-q (Figure lb). 

In addition to the conserved miRNAs, we also identified 
137 non-conserved miRNAs that belong to 44 miRNA 



families (Additional file 1) after removing miRNAs with 
expression levels that were too low to be analysed for dif- 
ferential expression. Among the 137non-conserved miR- 
NAs, some were legume-specific. The most abundant 
miRNA was miR2086 followed by miR2111 in the control 
(CK) library. In the drought stress (DS) library, miR2111 
was the most abundant, followed by miR1507. In addition, 
miR2630 was found to have the most members in the two 
libraries (miR2630a-y). 
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Figure 1 Distribution of small RNAs from control and drought stress libraries. The size distribution of small RNAs is shown in panel (a) 
with frequency. The absolute miRNA number in different conserved families is shown in panel (b). 
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Identification of new mIRNAs/new members of known 
miRNA families 

The formation of stable hairpin structure is one of the 
essential features for identification of new miRNAs [56]. 
To identify new miRNAs and new members of known 
miRNA families, we pooled the reads obtained from the 
two libraries, and identified the miRNAs by folding the 
sequences of potential miRNA precursors using the 
mfold web server [57]. We obtained 238 potential 
miRNA candidates (Table 2, Additional file 2). In addi- 
tion to the hairpin structure, detection of miRNA*s is 
another proof that has been widely used to confirm the 
new miRNAs [58]. Among the 238 potential miRNA 
candidates, we found 29 miRNAs with complementary 



miRNA*s (Table 2, Figure 2, Additional file 3), indicat- 
ing that these candidate miRNAs are likely to be new 
miRNAs or new members of known miRNA families in 
M. truncatula. The counts of several miRNA*s of these 
miRNAs were low (Table 2). This may make the anno- 
tation of these new miRNAs questionable. However, no 
requirement on the counts of several miRNA's was 
given in the most recent criteria for annotation of plant 
miRNAs [58]. Therefore, these miRNAs were still con- 
sidered to be new miRNAs or new members of known 
miRNAs. A similar identification of new miRNAs has 
also been reported by others in the literature. For 
instance, candidate miRNAs with only one miRNA* 
have been identified as new miRNAs in rice [59]. 



Table 2 New miRNAs and new members of known miRNA families with miRNA^s in M. truncatula 



mIRNAs 


Sequence (5'-3') 


Length 


Counts of miRNAs/ 


Location 


Arm 


Length of 


WIFE (kca 






(nt) 


miRNA*s 






precursors (nt) 


mol) 


miR5213 


uacgugugucuucaccucugaa 


22 


1 5045/44 


MtChr2:1 6737666:1 6737796: 
+ 


5' 


131 


-43.1 


miR5274b 


auaugacggaguguaaaugcc 


21 


29/2 


MtChr4:11424365:l 1424459:- 


5' 


95 


-663 


miR5554a 


ugugcaucuugaacaaugguau 


22 


134/3 


MtChrl :3 1 285750:3 1 285854:- 


5' 


105 


-52.8 


miR5554b 


ugugcaucuugaacaaugguau 


22 


134/3 


MtChr441 01 339341 01 3497:- 


5' 


105 


-573 


miR5554c 


ugugcaucuugaacaaugguau 


22 


134/1 


MtChr7:1 9742691 :1 9742795:- 


5' 


105 


-533 


miR5555 


uaagaguauaauaugacuuug 


21 


40664 


MtChrl :1 2905862:1 2905956:- 


5' 


95 


-47.5 


miR5556 


ugaugacqqaaqaaauccaaa 


21 


155/1 


MtChr4:231 99944:23200064:- 


3' 


121 


-54.8 


miR5557 


ugcuuccuuaguacuuguuga 


21 


5/3 


MtChr5:60261 67:6026481 :+ 


3' 


315 


-96.7 


miR5558 


uuuuccaauucuaagucuauc 


21 


139/8 


MtChr5:20521 087:20521 172: 


5' 


86 


-33 


miR5559 


uacuuggugaauuguuggauc 


21 


2124/2 


MtChr6:937306:937443:+ 


5' 


138 


-48.8 


miR5560 


ugccggcucaaugaaugcggag 


22 


62/2 


MtChr6:47474284747538:-F 


3' 


111 


-53.8 


miR5561 


cauuuggagagacauagacaa 


21 


449/3 


MtChr7:301 70572:301 70660:- 


5' 


89 


-42.9 


miR5562 


uguggagucuuuugcaugaag 


21 


16/1 


MtChr7:321 94444:321 94667:- 


5' 


224 


-564 


miR5563 


ugauaucaggcaacucggucc 


21 


12/1 


MtChr8:32699524:32699624:- 


5' 


101 


-52.9 


miR156j 


uugacagaagagggugagcaca 


22 


10341/23 


MtChrl:3228429:3228554:-F 


5' 


126 


-444 


miR167b 


ugaagcugccagcaugaucug 


21 


40937/1 


MtChr7:34009590:34009796: 


5' 


207 


-78.7 


miR168c 


ucgcuuggugcaggucgggaa 


21 


56024/1434 


MtChr5:1 0397699:1 0397824: 
4- 


5' 


126 


-68.8 


miR172b 


agaaucuugaugaugcugcau 


21 


35344/63 


AC235487:1 97944:1 98077:+ 


3' 


134 


-55.5 


miR172c 


agaaucuugaugaugcugcau 


21 


35386/8 


AC233663:1 01 68:1 0308:+ 


3' 


141 


-52.8 


miR408 


augcacugccucuucccuggc 


21 


75/55 


MtChrl21952074:21952198: 
+ 


3' 


125 


-45.1 


miR2mu 


uaaucugcauccugagguuua 


21 


413/23 


MtChr7:14331003:14331117:- 


5' 


115 


-50.2 


miR2111v 


uaaucugcauccugagguuua 


21 


1192/41 


MtChr7:1 4356392:1 4356499:- 


5' 


108 


-484 


miR2592a 


gaaaaacaugaaugucgagcg 


21 


28/1 


MtChrl 27380857:27381 050: 
+ 


3' 


194 


-894 


miR2592bl 


uggcaaguuugaauuuaccuca 


22 


43/1 


MtChr4:1 8239825:1 8239954: 
+ 


5' 


130 


-69.3 


miR2592bm 


ggaaaacaugaaugucgggug 


21 


918/294 


MtChr53041 01 2:3041 204:+ 


3' 


193 


-1044 


miR2592bn 


ggaaaacaugaaugucgggug 


21 


530/160 


MtChr5:84001 82:8400375:+ 


3' 


194 


-111.8 


miR2619b 


auauguuuugauucuuuggca 


21 


9/3 


MtChr4:6335670:6335840:+ 


5' 


171 


-94.7 


miR2543b 


uuugggaucagaaauuagaga 


21 


361/11 


MtChr541 83644941 836558:- 


3' 


110 


-38.8 


miR4414a 


agcugcugacucguugguuca 


21 


663/45 


MtChrl 305 18803:305 18924: 


5' 


122 


-50.5 



+ 
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miR5555 

10 20 30 40 

UAC A — -I UU 

AUGUAUCUAAGAGU UAAUAUGACUUUGAUUAUG ACA CAC \ 
UACAUAGAUUCUCA AUUAUGCUGAAGCUAAUAC UGU GUG U 
AAC C CGU A" AG 

90 80 70 60 50 



miR5558 

10 20 30 40 

U| CAUA C CCUAU C 

GACCU UUUUC AAUUCUAAGUCUAUC GAAUUC A 
CUGGA AAAAG UUAAGAUUUAGAUAG CUUAAG A 
C" AUCG A AGC~ U 

80 70 60 50 



miR5561 

10 20 30 40 

-I A C AA U G 

UGGA ACUUUCAUUUGGAGAGA AUAGACA UGAAAU AAAU A 
ACCU UGAAGGUAAAUCUCUCU UAUCUGU ACUUUA UUUA U 
U" G C A- - A 

80 70 60 50 



miR5563 

10 20 30 40 

UUA C AC — I CA 

CAAUU AAUGAUAU AGGCAACUCGGUCCUUCUGUUA UUGGC \ 
GUUAA UUGCUGUA UCCGUUGAGUCAGGAAGGCAAU AACUG A 
G UAC A GACCA" UC 

90 80 70 60 50 



miR167b 

10 20 30 40 50 60 70 80 90 100 110 

CAAl U G G U UUCCUGUGU U UAAU CUGA C A A U UAA 

UGACAG UGAAGCU CCAGCAUGAUCUGU CUU CC AUA ACUU UUCCAG AUUUAAAUAUAA CAAAAA AUA AUA GUUUGGUC A 

AUUGUC ACUUCGA GGUUGUACUAGACA GAA GG UAU UGAG GGGGUC UGAAUUUGUAUU GUUUUU UAU UAU CAAACUAG U 

UAA" C - A U UGUU U U — C A A U UUU 

200 190 180 170 160 150 140 130 120 

Figure 2 Representatives of predicted precursors' hairpin structures of new miRNA/new members of known miRNA families The 
mature miRNA and miRNA* sequences are coloured in red and blue, respectively. All the predicted hairpin structures of these mlRNA precursors 
are shown in Additional file 3. 



Twenty-two out of 29 new miRNAs/new members of 
known miRNA families had length of 21 nt, while 
remaining miRNAs had length of 22 nt. The majority of 
the miRNAs started with a 5' uridine, a hallmark of miR- 
NAs [60]. The minimum free energy (MFE) for hairpin 
structure of miRNA precursors was lower than -30 kcal/ 
mol (Table 2). This feature is in agreement with other 
reported values in the literature [61]. The length of new 
miRNAs' precursors ranged from 86 nt to 315 nt (Table 
2). These values are comparable to those reported for M. 
truncatula in the literature [42,51,52,62,63]. 

Response of known miRNAs to drought stress 

To identify drought-responsive miRNAs, the normalized 
expression of miRNAs in the two libraries (control, CK 
and drought stress, DS) was compared. Based on the 
results of high-throughput sequencing, we selected those 
miRNAs with changes in expression levels being greater 
than 1.5-fold in response to drought treatment (Figure 3a 
and 3b, Additional file 4) to validate the expression pat- 
terns by real-time quantitative PCR. As shown in Figure 3b 
and 3c, expression patterns of drought-responsive miRNAs 
from high-throughput sequencing and RT-qPCR exhibited 
similar results. Results from both the methods showed that 
22 members in 4 miRNA families, i.e., miR399, miR2089, 
miR2111 and miR2118, were up-regulated in response to 



the drought stress. Conversely, 10 members belonging to 6 
miRNA families, i.e., miR164, miR169, miR171, miR396, 
miR398 and miRlSlO, were down-regulated in response to 
drought stress (Figure 3b and 3c). These miRNAs have 
been reported to be involved in diverse cellular processes 
in plants [14,17,22,52]. The known target genes of these 
miRNAs and their function annotations were summarized 
in Table 3. For those miRNAs whose targets are not 
known, we predicted their targets using the psRNATarget 
http://bioinfo3.noble.org/psRNATarget/ and the srna- tools 
http://srna-tools.cmp.uea.ac.uk/plant/ [64], and the results 
are given in the Table 4. Among these miRNAs, several 
miRNAs have been reported to be involved in abiotic stres- 
ses. For example, miR399 and miR2111 have been reported 
to be up-regulated by phosphate starvation [23,65], while 
miR169 with target of CCAAT Binding Factor is down- 
regulated in response to drought stress [14,17]. 

Response of new miRNAs/new members in known miRNA 
families to drought stress 

We further examined whether the 29 new miRNAs/new 
members of known miRNA families were responsive to 
drought stress using the same method as the known 
miRNAs. The miRNAs with changes in expression levels 
being greater than 1.5-fold and /7-values less than 0.05 
are presented in Table 4. Results from high-throughput 
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Normalized expression level of control miRNAs 



Figure 3 Differential expression analysis between control and drought stress. Data points at tine upper and lower side of tine slope line 
represent up- and down-regulated miRNAs in panel (a), respectively. The changes in miRNAs for up- and down-regulated ones are greater than 
1.5-fold. Other miRNAs include those that are not responsive to drought stress and those changes induced by drought stress at p > 0.05. 
Expression of control and drought stress was normalized on the basis of 1 M reads. Differential expression of known miRNAs in response to 
drought stress is shown in panel (b). The positive and negative values mean miRNAs whose expression was stimulated and suppressed by 
drought stress, respectively. * and ** mean significant difference between control and drought stress at 0.01 < p < 0.05 and p < 0.01, 
respectively. The relative expression level of miRNAs measured by RT-qPCR in response to drought stress is shown in panel (c). 



sequencing and RT-qPCR showed that, among the 14 
new miRNAs, 5 miRNAs were responsive to drought 
stress with miR5274b and miR5558 being up-regulated, 
while miR5554a-c being down-regulated. In addition, 
among the 15 new members of known miRNA families, 
3 miRNAs were responsive to drought stress with 
miR2111u, V being up-regulated and miR4414a being 
down-regulated (Figure 3b and 3c, Table 4). 

Discussion 

High-throughput sequencing of M. truncatula small RNAs 

High-throughput sequencing has been used to study miR- 
NAs at whole genome level in several plant species. 



including Arabidopsis [39,40], rice [43], wheat [44], 
soybean [45], cotton [46], grapevine [53], Medicago 
[42,51,52], Brachypodium [20], tomato [41], populus [47], 
citrus [48], peanut [49], Porphyra [50]. However, most of 
the studies on high-throughput sequencing miRNAs have 
focused on miRNAs under non-stressed, normal growth 
conditions, and only a few studies have compared miRNAs 
under control conditions to those under conditions of 
abiotic stresses by high-throughput sequencing plant miR- 
NAs [20,43,44]. For instance, Sunkar et al. (2008) reported 
a total of 714,202 reads in rice seedlings from three inde- 
pendent libraries (control, drought and salt stress libraries) 
[43]. Twenty-eight and 12 miRNAs in response to cold 



Table 3 The known targets of drought-responsive miRNAs and their function annotations 


miRNA 
families 


Expression 
pattern 


Targets 


Functions and responsiveness 


References 


miR164 


down 


NAC domain transcription 
factors 


lateral root development 


[9] 


miR159 


down 


CCAAT Binding Factor (CBF) 


response to drought, cold and salinity, nodule development 


[14,15,17,18,22,33,74] 


miR171 


down 


GRAS transcription factors 


response to drought, cold and salinity, nodule 
morphogenesis, floral development 


[5,18,33] 


miR396 


down 


Growth Regulating Factor (GRF) 


response to drought and salinity, cell proliferation 


[18,52,75] 


miR398 


down 


Cu/Zn superoxide dismutases 
(CSD1, CSD2) 


response to oxidative stress 


[15,27] 


miR399 


up 


PH02/ubiquitin conjugating 
enzyme 


balance of phosphorus 


[23] 


miR2118 


up 


TIR-NBS-LRR domain protein 


response to drought cold, salinity, and ABA 


[51,76] 
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Table 4 Predicted targets of drought-responsive miRNAs 



miRNAs 



Expression 
pattern 



Predicted targets 



miRISlOa 



miR2089 

miR2111a- 
s 

miR21 1 1u, 



miR5274b 

miR5554a- 

c 

miR5558 



down 1. pyruvate decarboxylase (PDC) isozyme 

1 

2. concanavalin A-like lectin/glucanase 
3. F-box protein 
up NB-ARC domain protein 

up 1. calcineurin-like phosplioesterase 

2. membrane protein SAK 



up 



up 
down 

up 



1. calcineurin-like phosphoesterase 
family protein 

2. membrane protein SAK 

DNA-damage-repair/toleration protein 

polynucleotidyl transferase, Ribonuclease 
H fold 

1. initiation factor elF-4 gamma 

2. homeodomain-related POX 



Stress in Brachypodium [20,44] and heat stress in wheat 
[44] have been identified by high-throughput sequencing, 
respectively. In the present study, we constructed two 
sRNA hbraries, non-stressed, control and drought stressed 
in legume model plant, M. truncatula and identified 
known and new miRNAs that were responsive to drought 
stress. 

The successful application of high-throughput sequen- 
cing technology to systemically identify plant miRNAs has 
greatly advanced our knowledge on the functions of miR- 
NAs in plants in recent years. There have been several 
reports on the identification of miRNAs in M. truncatula 
in the literature. For instance, 26,656 and 844,110 sRNA 
reads were reported in two recent studies from M. trunca- 
tula by Jagadeeswaran et al. (2009) and Lelandais-Briere et 
al. (2009), respectively, with 97,028 unique sequences 
being obtained by the latter studies [51,52]. However, in 
both cases, sRNAs were sequenced using Roche 454 
sequencer, and the total reads and unique sequences 
obtained from these studies are much less than we 
reported in the present study (Table 1). In a study by Szit- 
tya et al. (2008), the authors obtained 3,948,871 reads and 
1,563,959 unique sequences from two sRNAs libraries of 
M. truncatula using the same high-throughput sequencing 
technology (Illumina-Solexa) [42]. In contrast to these 
reported findings, we got a total of 13,683,619 reads and 
4,878,445 unique sequences from two libraries of M. trun- 
catula. Therefore, our results contain more reads than 
those previously reported for miRNAs in M. truncatula in 
the literature. Moreover, the high-throughput sequencing 
technology used in the present study allows us to identify 
the miRNAs with low abundance, thus accounting for the 
greater unique sequences found in our studies. In addition 
to the greater amounts of miRNAs sequencing data, the 



use of most updated miRBase also contributes to our iden- 
tification of more known miRNAs in M. truncatula. The 
available database for M. truncatula contains more miR- 
NAs than those for other plant species. For example, in 
the database, there are 19 miRNAs and 41 miRNAs for 
Brachypodium and wheat, respectively, while there are 375 
miRNAs in the database for M. truncatula (miRBase 17). 
In addition, to the best of our knowledge, we are the first 
one to use the M. truncatula genome sequence Mt3.0 to 
analyze miRNAs in this species. The more genomic infor- 
mation in Mt3.0 than Mt2.0 may also contribute to pre- 
dicting more new miRNAs in M. truncatula. 

Drought-responsive miRNAs 

Plant miRNAs have a strong propensity to target genes 
associated with development, particularly those genes 
encoding transcription factors and F-box proteins [25] . In 
the present study, we found that the known or predicted 
targets of miR164, miR169, miR171, miR396, miR1510 
and miR5558 were either transcription factors or F-box 
proteins (Table 3, 4). Under drought stress, an increase in 
root/shoot ratio was found in M. truncatula (data not 
shown). miR164 has been reported to regulate root devel- 
opment by a homeostatic mechanism to clear NACl 
mRNA, leading to down-regulation of auxin signals. It has 
been shown that decreases in NACl mRNA levels due to 
inducible expression of miR164, thus resulting in reduc- 
tion in lateral root emergence in Arabidopsis [9]. There- 
fore, it is conceivable that the suppression of miR164 
expression may contribute to the increase in root/shoot 
ratio under drought stress. 

Several studies have revealed that miR169 is responsive 
to abiotic stresses such as drought, cold and salinity in 
different species [14,15,17,22,33]. Under drought stress, 
miR169 exhibited different expression patterns among dif- 
ferent species. For example, in Arabidopsis, miR169 is 
down-regulated by drought stress through an ABA-depen- 
dent pathway, resulting in accumulation of NFYA5 with 
high affinity and sequence specificity for the CCAAT box, 
which is crucial for the expression of a number of 
drought-responsive genes [14]. In contrast to Arabidopsis, 
miR169g in rice is up-regulated by drought [17,18]. In 
rice, DREs (dehydration-responsive element) as the 
upstream of MIR169g are regulated by cold and drought 
stress, thus leading to up-regulation of miR169g [17]. 
Trindade et al. (2010) reported that the expression of 
miR169 in leaves of M. truncatula is not responsive to 
drought stress [15]. In contrast, our results from both the 
high-throughput sequencing and RT-qPCR showed that 
miR169 was down-regulated under drought stress in 
M. truncatula (Figure 3b and 3c, Additional file 4). This 
discrepancy may result from the difference in degree of 
drought stress exposed to plant materials in the two stu- 
dies. For instance, in the present study, the drought 
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Stressed samples were collected after exposing of plants to 
drought for varying period (6, 8, 10 and 12 d after with- 
holding water), thus our samples included plant materials 
suffering from a wide range of drought stress. This type of 
stress has been used in study of molecular response of 
plants to drought stress [66,67] . In the studies of Trindade 
et al. (2010), water stress was exposed by suppressing 
water supply with relative water content in leaves being 
approx. 50% and 30%, respectively [15]. In the present 
paper, the drought stressed leaves with relative water con- 
tent of 87.4%, 78.8%, 75.2% and 68.3% were used for analy- 
sis of miRNAs. Therefore, our samples were mainly those 
suffering from mild and moderate water stress compared 
to those of Trindale et al. (2010). If miR169 is responsive 
to mild and/or early drought stress exclusively, the respon- 
siveness of miR169 may not be detected by the more 
severe drought stress used by the authors [15]. In addition 
to withholding water supply, several studies on effect of 
drought on miRNAs also treated plants with PEG or man- 
nitol for varying period [17,33]. It is difficult to compare 
the natural drought stress with PEG-induced and/or man- 
nitol osmotic stress as the two treatments may differ in 
induction of water stress in terms of rapidity and severity. 

Accumulation of reactive oxygen species (ROS) is a 
common phenomenon in response of plants to abiotic 
stress. The accumulated ROS damage nucleic acid, oxidize 
proteins and cause lipid peroxidation [68,69]. Superoxide 
dismutases (SODs) detoxify superoxide radicals. The tar- 
gets of miR398 are two Cu/Zn superoxide dismutases 
(cytosolic CSDl and chloroplastic CSD2), and miR398 
expression was reported to be down-regulated transcrip- 
tionally by oxidative stresses [27] . Oxidative stress often 
occurs concurrently with drought stress. In the present 
study, we found that miR398 was down-regulated under 
drought stress. This would lead to increases in activities of 
SODs. The drought-induced down-regulation of miR398 
in M. truncatula is consistent with the results in maize 
[16], but it is contrast to the results reported by Trindale 
etal. (2010) [15] and Kantar et al. (2011) [34]. The differ- 
ences in the expression of miR398 among different studies 
may results from differences in species, extent and dura- 
tion of drought stress in different studies. 

Plants suffering from water-deficit often display reduced 
uptake of mineral nutrients. In this context, miR399 nega- 
tively regulates the concentration of inorganic phosphate 
(Pi) by targeting PH02, a type of E2 conjugase, and over- 
expression of miR399 in Pi-replete conditions represses 
expression of E2 conjugase, leading to an increase in Pi 
concentration in leaves in Arabidopsis [23] . Expression of 
miR399 is reduced under Pi-deprived conditions to facili- 
tate accumulation of Pi in plants. It has been verified that 
miR2111 is up-regulated by Pi starvation [65]. In our 
study, the expression miR399 and miR2111 was similar 
under drought stress. These results indicate that the 



function of miR2111 may be as important as miR399 in 
regulation of nutrient acquisition. 

Exposure of plants to a moderate stress induces resis- 
tance to other stresses, a phenomenon known as cross 
adaptation, which has been found in different combina- 
tions of stresses [70]. For example, osa-miR821 isolated 
from virus-infected rice tissues is also expressed in roots 
of salt-stressed plants, while it is not expressed in healthy, 
non-stressed plants [71]. Under drought stress, we found 
that there was up-regulation of miR2089 and miR2118, 
whose targets may be proteins associated with disease 
resistance. It is envisaged that these miRNAs may enhance 
the ability of drought tolerance through unknown 
mechanisms associated with cross adaptation in plants. 
Future work aiming at functional elucidation of these miR- 
NAs is warranted by over-expressing these miRNAs in 
M. truncatula. 

Conclusions 

We obtained a total of 13,683,619 reads from two small 
RNA libraries of M. truncatula by high-throughput 
sequencing. Twenty-two members in 4 miRNA families 
and 10 members belonging to 6 miRNA families were 
found to be up-regulated and down-regulated in response 
to drought stress by both high-throughput sequencing and 
RT-qPCR. In addition, we also predicted 29 new miRNAs/ 
new members of known miRNA families, of which 8 miR- 
NAs were responsive to drought stress by high-throughput 
sequencing and RT-qPCR. These findings provide valuable 
information for further functional characterization of miR- 
NAs in response to abiotic stress in general and drought 
stress in particular in legume plants. 

Methods 

Plant materials and drought stress treatment 

Four Medicago truncatula (cv Jemalong A17) seedlings 
were grown in a pot (diameter 10 cm) filled with vermicu- 
lite: peat soil (2:1) under controlled conditions (26°C day/ 
20°C night, 14-h photoperiod, and 50% relative humidity). 
Drought stress was initiated by withholding water supply 
to 3-week-old seedlings for varying periods after seedlings 
were fully watered. The severely wilted leaves appeared on 
the 14th day after the water withholding. The soil water 
content was reduced from approx. 60% to 8.1% during the 
drought stress. Shoots under drought stress were har- 
vested after withholding of water for 6, 8, 10 and 12 d, and 
mbced averagely as drought treatment materials (DS). The 
relative water content for leaves collected after 6, 8, 10 and 
12 d drought stress was 87.4%, 78.8%, 75.2% and 68.3%, 
respectively. At the same time, shoots of M. truncatula 
seedlings grown under normal watering conditions were 
also harvested as control (CK). The relative water content 
for control leaves was 91.4%. The measure method of rela- 
tive water content was described by Catsky (1960) [72]. 
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Small RNA libraries construction for high-throughput 
sequencing 

To construct small RNA libraries, total RNA was extracted 
from shoots of control (CK) and drought stress (DS) using 
the Trizol (Invitrogen) according to the manufacturer's 
protocols. For each sample, the 18-30 nt small RNAs were 
ligated with 5' and 3' RNA adapter by T4 RNA ligase 
(TaKaRa) after they were purified by electrophoretic 
separation on a 15% TBE-urea denaturing PAGE gel, and 
at each step purified by urea PAGE gel electrophoretic 
separation. The RNA was subsequently transcribed to sin- 
gle-stranded cDNA using Superscript II reverse transcrip- 
tase (Invitrogen). Thereafter the cDNA was used as 
templates for double-stranded cDNA synthesis by PGR 
amplification using primers that anneal to adapters. The 
purified DNA was sequenced on a Solexa sequencer (Illu- 
mina). The raw data have been submitted to Gene Expres- 
sion Omnibus (GEO, http://www.ncbi.nlm.nih.gov/ 
projects/geo/) and the accession number is GSE29154. 



sample (Normalized expression = Actual miRNA count/ 
Total count of clean reads*l,000,000). miRNAs whose 
normalized expression of two libraries is smaller than 
one were removed, because their expression levels are 
too low. The fold change between treatment and control 
was calculated as: Fold-change = log2 (DS/CK). The sta- 
tistical analysis was performed according to Poisson dis- 
tribution. The p-value was calculated by the following 
formula. 



{x + y) ! 



No 

p(x\y) = i—Y- , 

^ ^ Ni 
C(y < yminix) = ^ p(y|x) 

oo 

D(y > ymaxlx) = ^ p(y|x) 



Bioinformatics analysis 

After the sequence tags from Solexa sequencing went 
through the data cleaning by removing the low quality 
tags (i.e., tags less than 18 nt and tags whose adaptors 
were null) and contaminants (adaptors and polyA), the 
length distribution and common/specific sequences 
between two samples were analyzed. After removing adap- 
tor sequences of the left high-quality small RNA reads 
with exact matches to the adaptor sequences, reads were 
mapped to the M. truncatula genome sequence (MtS.O) 
downloaded from the website http://www.medicagohap- 
map.org/downloads_genome/Mt3 using SOAP [54]. 
rRNAs, tRNAs, snRNAs and snoRNAs were removed 
from the matched sequences through BLASTn search [73] 
using NCBI Genebank database http://www.ncbi.nlm.nih. 
gov/blast/Blast.cgi/ and Rfam database http://www.sanger. 
ac.uk/Software/Rfam/. Mismatches were not allowed in 
the above two approaches. The unique sequences left were 
aligned with known miRNAs from miRBase 17 http:// 
www.mirbase.org/ [55]. The potential candidate miRNAs 
were identified by folding the flanking genome sequence 
of unique small RNAs using the mfold web server [57] . 
Parameters were set based on the criteria for annotation 
plant miRNAs [58]. Target predictions were performed 
using the psRNATarget http://bioinfo3.noble.org/psRNA- 
Target/ and the srna-tools http://srna-tools.cmp.uea.ac.uk/ 
plant/ [64], through aligned with genome of M. truncatula 
and A. thaliana. 



Real-time quantitative PCR of mature miRNAs 

RT-qPCR was used to validate the results obtained from 
the high-throughput sequencing of miRNAs. RNA that 
was isolated using the Trizol (Invitrogen) as described 
above was reversely transcribed using One Step Prime- 
Script miRNA cDNA Synthesis Kit (TaKaRa). This kit 
adds ploy (A) to the 3' end of miRNAs and start to 
reversely transcribe. The reverse transcription was led 
by a kind of special oligo-dT ligated with known 
sequence at its 5' end. RT-qPCR was performed using 
SYBR Premix Ex Tag II (TaKaRa) and all the primers 
used were listed in Additional file 5. Small nuclear RNA 
U6 was used as internal control. Subsequently, RT- 
qPCR was performed using Stratagene M x 3000P 
instrument. 

Additional material 



Additional file 1: The miRNA abundance of non-conserved miRNA 
families in control (CK) and drought stress (DS) libraries. 

Additional file 2: Potential miRNA candidates without miRNA*s In 
M. truncatula 

Additional file 3: The predicted hairpin structures of all the 29 new 
miRNAs/new members of known miRNA families' precursors 

Additional file 4: Known miRNAs in response to drought stress 

Additional file 5: The primers designed for RT-qPCR 



Differential expression analysis of miRNAs under the 
drought stress 

The frequency of miRNAs of two libraries was normal- 
ized to one million by total number of miRNAs in each 
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